Example use case: Common-envelope evolution

In this notebook we look at how common-envelope evolution (CEE) alters binary-star orbits. We construct a population of low- and intermediate-mass binaries and compare their orbital periods before and after CEE. Not all stars evolve into this phase, so we have to run a whole population to find those that do. We then have to construct the pre- and post-CEE distributions and plot them.

First, we import a few required Python modules.

[1]:
import os
import math
import matplotlib.pyplot as plt
from binarycpython.utils.functions import temp_dir, output_lines
from binarycpython import Population

TMP_DIR = temp_dir("notebooks", "notebook_comenv", clean_path=True)
VERBOSITY = 2

Setting up the Population object

We set up a new population object. Our stars evolve to \(13.7\text{ }\mathrm{Gyr}\), the age of the Universe, and we assume the metallicity \(Z=0.02\). We also set the common-envelope ejection efficiency \(\alpha_\mathrm{CE}=1\) and the envelope structure parameter \(\lambda=0.5\). More complex options are available in binary_c, such as \(\lambda\) based on stellar mass, but this is just a demonstration example so let’s keep things simple.

[2]:
# Create population object
population = Population(tmp_dir=TMP_DIR, verbosity=VERBOSITY)
population.set(
    # grid options
    log_dt = 10, # log progress every 10 seconds

    # binary-star evolution options
    max_evolution_time=13700,  # maximum stellar evolution time in Myr (13700 Myr == 13.7 Gyr)
    metallicity=0.02, # 0.02 is approximately Solar metallicity
    alpha_ce = 1.0,
    lambda_ce = 0.5,
)
adding: log_dt=10 to population_options
adding: max_evolution_time=13700 to BSE_options
adding: metallicity=0.02 to BSE_options
adding: alpha_ce=1.0 to BSE_options
adding: lambda_ce=0.5 to BSE_options

Stellar Grid

We now construct a grid of stars, varying the mass from \(1\) to \(6\text{ }\mathrm{M}_\odot\). We avoid massive stars for now, and focus on the (more common) low- and intermediate-mass stars. We also limit the period range to \(10^4\text{ }\mathrm{d}\) because systems with longer orbital periods will probably not undergo Roche-lobe overflow and hence common-envelope evolution is impossible.

[3]:
# Set resolution and mass range that we simulate
resolution = {"M_1": 10, "q" : 10, "per": 10}
massrange = [1, 6]
logperrange = [0.15, 4]

population.add_sampling_variable(
    name="lnm1",
    longname="Primary mass",
    valuerange=massrange,
    samplerfunc="self.const_linear(math.log({min}), math.log({max}), {res})".format(min=massrange[0],max=massrange[1],res=resolution["M_1"]),
    precode="M_1=math.exp(lnm1)",
    probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
    dphasevol="dlnm1",
    parameter_name="M_1",
    condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
)

# Mass ratio
population.add_sampling_variable(
     name="q",
     longname="Mass ratio",
     valuerange=["0.1/M_1", 1],
     samplerfunc="self.const_linear({}/M_1, 1, {})".format(massrange[0],resolution['q']),
     probdist="self.flatsections(q, [{{'min': {}/M_1, 'max': 1.0, 'height': 1}}])".format(massrange[0]),
     dphasevol="dq",
     precode="M_2 = q * M_1",
     parameter_name="M_2",
     condition="",  # Impose a condition on this grid variable. Mostly for a check for yourself
)

# Orbital period
population.add_sampling_variable(
    name="log10per", # in days
    longname="log10(Orbital_Period)",
    valuerange=[0.15, 5.5],
    samplerfunc="self.const_linear({}, {}, {})".format(logperrange[0],logperrange[1],resolution["per"]),
    precode="""orbital_period = 10.0 ** log10per
sep = calc_sep_from_period(M_1, M_2, orbital_period)
sep_min = calc_sep_from_period(M_1, M_2, 10**{})
sep_max = calc_sep_from_period(M_1, M_2, 10**{})""".format(logperrange[0],logperrange[1]),
    probdist="self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**{}), math.log10(10**{}), {})".format(logperrange[0],logperrange[1],-0.55),
    parameter_name="orbital_period",
    dphasevol="dlog10per",
)
Added sampling variable: {
    "name": "lnm1",
    "parameter_name": "M_1",
    "longname": "Primary mass",
    "valuerange": [
        1,
        6
    ],
    "samplerfunc": "self.const_linear(math.log(1), math.log(6), 10)",
    "precode": "M_1=math.exp(lnm1)",
    "postcode": null,
    "probdist": "self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
    "dphasevol": "dlnm1",
    "condition": "",
    "gridtype": "centred",
    "branchpoint": 0,
    "branchcode": null,
    "topcode": null,
    "bottomcode": null,
    "sampling_variable_number": 0,
    "dry_parallel": false,
    "dependency_variables": null
}
Added sampling variable: {
    "name": "q",
    "parameter_name": "M_2",
    "longname": "Mass ratio",
    "valuerange": [
        "0.1/M_1",
        1
    ],
    "samplerfunc": "self.const_linear(1/M_1, 1, 10)",
    "precode": "M_2 = q * M_1",
    "postcode": null,
    "probdist": "self.flatsections(q, [{'min': 1/M_1, 'max': 1.0, 'height': 1}])",
    "dphasevol": "dq",
    "condition": "",
    "gridtype": "centred",
    "branchpoint": 0,
    "branchcode": null,
    "topcode": null,
    "bottomcode": null,
    "sampling_variable_number": 1,
    "dry_parallel": false,
    "dependency_variables": null
}
Added sampling variable: {
    "name": "log10per",
    "parameter_name": "orbital_period",
    "longname": "log10(Orbital_Period)",
    "valuerange": [
        0.15,
        5.5
    ],
    "samplerfunc": "self.const_linear(0.15, 4, 10)",
    "precode": "orbital_period = 10.0 ** log10per\nsep = calc_sep_from_period(M_1, M_2, orbital_period)\nsep_min = calc_sep_from_period(M_1, M_2, 10**0.15)\nsep_max = calc_sep_from_period(M_1, M_2, 10**4)",
    "postcode": null,
    "probdist": "self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**4), -0.55)",
    "dphasevol": "dlog10per",
    "condition": null,
    "gridtype": "centred",
    "branchpoint": 0,
    "branchcode": null,
    "topcode": null,
    "bottomcode": null,
    "sampling_variable_number": 2,
    "dry_parallel": false,
    "dependency_variables": null
}

Logging and handling the output

We now construct the pre- and post-common envelope evolution data for the first common envelope that forms in each binary. We look at the comenv_count variable, we can see that when it increases from 0 to 1 we have found our object. If this happens, we stop evolution of the system to save CPU time.

[4]:
custom_logging_statement = """

/*
 * Detect when the comenv_count increased
 */
if(stardata->model.comenv_count == 1 &&
   stardata->previous_stardata->model.comenv_count == 0)
{
   /*
    * We just had this system's first common envelope:
    * output the time at which this happens,
    * the system's probability (proportional to the number of stars),
    * the previous timestep's (pre-comenv) orbital period (days) and
    * the current timestep (post-comenv) orbital period (days)
    */
    Printf("COMENV %g %g %g %g\\n",
           stardata->model.time,
           stardata->model.probability,
           stardata->previous_stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS,
           stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS);

    /*
     * We should waste no more CPU time on this system now we have the
     * data we want.
     */
    stardata->model.evolution_stop = TRUE;
}
"""

population.set(
    C_logging_code=custom_logging_statement
)
adding: C_logging_code=

/*
 * Detect when the comenv_count increased
 */
if(stardata->model.comenv_count == 1 &&
   stardata->previous_stardata->model.comenv_count == 0)
{
   /*
    * We just had this system's first common envelope:
    * output the time at which this happens,
    * the system's probability (proportional to the number of stars),
    * the previous timestep's (pre-comenv) orbital period (days) and
    * the current timestep (post-comenv) orbital period (days)
    */
    Printf("COMENV %g %g %g %g\n",
           stardata->model.time,
           stardata->model.probability,
           stardata->previous_stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS,
           stardata->common.orbit.period * YEAR_LENGTH_IN_DAYS);

    /*
     * We should waste no more CPU time on this system now we have the
     * data we want.
     */
    stardata->model.evolution_stop = TRUE;
}
 to population_options

The parse function must now catch lines that start with “COMENV” and process the associated data. We set up the parse_data function to do just this.

[5]:
from binarycpython.utils.functions import bin_data, datalinedict
import re

# log-period distribution bin width (dex)
binwidth = 0.5

def parse_function(self, output):
    """
    Parsing function to convert HRD data into something that Python can use
    """

    # list of the data items
    parameters = ["header", "time", "probability", "pre_comenv_period", "post_comenv_period"]

    # Loop over the output.
    for line in output_lines(output):
        # obtain the line of data in dictionary form
        linedata = datalinedict(line,parameters)

        # choose COMENV lines of output
        if linedata["header"] == "COMENV":
            # bin the pre- and post-comenv log10-orbital-periods to nearest 0.5dex
            binned_pre_period = bin_data(math.log10(linedata["pre_comenv_period"]), binwidth)

            # but check if the post-comenv period is finite and positive: if
            # not, the system has merged and we give it an aritifical period
            # of 10^-100 days (which is very much unphysical)
            if linedata["post_comenv_period"] > 0.0:
                binned_post_period = bin_data(math.log10(linedata["post_comenv_period"]), binwidth)
            else:
                binned_post_period = bin_data(-100,binwidth) # merged!

            # make the "histograms"
            self.population_results['pre'][binned_pre_period] += linedata["probability"]
            self.population_results['post'][binned_post_period] += linedata["probability"]

    # verbose reporting
    #print("parse out results_dictionary=",self.population_results)

# Add the parsing function
population.set(
    parse_function=parse_function,
)
adding: parse_function=<function parse_function at 0x7fba7a927e50> to population_options

Evolving the grid

Now we actually run the population. This may take a little while. You can set num_cores higher if you have a powerful machine.

[6]:
# set number of threads
population.set(
    # set number of threads (i.e. number of CPU cores we use)
    num_cores=4,
    verbosity=VERBOSITY
)

# Evolve the population - this is the slow, number-crunching step
analytics = population.evolve()

# Show the results (debugging)
#print (population.population_results)
adding: num_cores=4 to population_options
adding: verbosity=2 to population_options
Creating the code for the shared library for the custom logging
File/directory /tmp/binary_c_python-david/notebooks/notebook_comenv/custom_logging/custom_logging.c doesn't exist. Can't remove it.
Writing the custom logging code to /tmp/binary_c_python-david/notebooks/notebook_comenv/custom_logging/custom_logging.c
File/directory /tmp/binary_c_python-david/notebooks/notebook_comenv/custom_logging/libcustom_logging_c2c13e8068bc4385b960647d0028913e.so doesn't exist. Can't remove it.
Calling the binary_c config code to get the info to build the shared library
Got options to compile:
        cc = cc
        ccflags = -std=gnu18 -DOPERATING_SYSTEM=linux -DLINUX -DPOSIX -DLARGEFILE_SOURCE -DMEMORY_ALLOCATION_MODEL=MEMORY_ALLOCATION_NATIVE -UALIGNSIZE -fstrict-aliasing -Wstrict-aliasing -g -Wno-sizeof-pointer-div -Wpedantic -Wshadow -Wno-variadic-macros -fstack-protector-all -rdynamic -fsignaling-nans -march=native -mtune=native -frounding-math -fno-stack-protector -ffloat-store -D__ACCURATE_BINARY_C__ -fno-finite-math-only -fasynchronous-unwind-tables  -export-dynamic -O0 -DCPUFREQ=4500 -DBINARY_C_SRC=/home/david/projects/binary_c_root/binary_c/src -DPGO_off -DBINUTILS_VERSION=2.34 -DBFD_2_33 -D_FILE_OFFSET_BITS=64 -D__HAVE_LINK_H -D__HAVE__VA_OPT__ -D__HAVE_GNU_QSORT_R -D__HAVE_NATIVE_EXP10 -D__HAVE_POSIX_FADVISE -DFPU_CONTROL -D__HAVE_ATTRIBUTE___RESTRICT____ -D__HAVE_ATTRIBUTE_ALLOC_SIZE__ -D__HAVE_ATTRIBUTE_AUTO_TYPE__ -D__HAVE_ATTRIBUTE_BUILTIN_EXPECT__ -D__HAVE_ATTRIBUTE_BUILTIN_ISNAN__ -D__HAVE_ATTRIBUTE_CONST__ -D__HAVE_ATTRIBUTE_DEPRECATED__ -D__HAVE_ATTRIBUTE_GNU_PRINTF__ -D__HAVE_ATTRIBUTE_HOT__ -D__HAVE_ATTRIBUTE_MALLOC__ -D__HAVE_ATTRIBUTE_NONNULL__ -D__HAVE_ATTRIBUTE_NORETURN__ -D__HAVE_ATTRIBUTE_PACKED__ -D__HAVE_ATTRIBUTE_PURE__ -D__HAVE_ATTRIBUTE_RETURNS_NONNULL__ -D__HAVE_ATTRIBUTE_UNUSED__ -DGIT_REVISION=6895:20230518:1fb08c5af -DGIT_URL=git@gitlab-dhendriks:binary_c/binary_c.git -DGIT_BRANCH=releases/2.2.4 -D__HAVE_LIBC__ -D__HAVE_LIBCFITSIO__ -D__HAVE_LIBGSL__ -I/home/david/gsl/include -D__HAVE_LIBGSLCBLAS__ -D__HAVE_LIBDL__ -D__HAVE_LIBPTHREAD__ -D__HAVE_LIBUUID__ -D__HAVE_LIBZ__ -D__HAVE_LIBBACKTRACE__ -D__HAVE_LIBBZ2__ -D__HAVE_LIBCDICT__ -D__HAVE_LIBM__ -D__HAVE_LIBRINTERPOLATE__ -D__HAVE_BACKTRACE_H__ -D__HAVE_IEEE754_H__ -D__HAVE_DRAND48__ -D__HAVE_HSEARCH_DATA__ -D__HAVE_MALLOC_H__ -D__HAVE_SETITIMER__ -D__HAVE_HAS_INCLUDE -D__HAVE_PKG_CONFIG__ -D__HAVE_VALGRIND__ -D__SHOW_STARDATA__ -D__DIFF_STARDATA__ -D__HAVE_7Z__ -D__HAVE_BZCAT__ -D__HAVE_ZCAT__ -O0 -shared -D_SEARCH_H
        ld = ld
        libs = -L/home/david/projects/binary_c_root/binary_c/src -L/usr/lib/x86_64-linux-gnu -L/home/david/gsl/lib -L/home/david/projects/binary_c_root/robs_extra_packages/lib -L/home/david/projects/binary_c_root/binary_c/src -lc -lcfitsio -lpthread -lgsl -lgslcblas -lm -ldl -luuid -lz -lbacktrace -lbz2 -lcdict -lrinterpolate -lbinary_c
        inc = -I/home/david/projects/binary_c_root/binary_c -I/home/david/projects/binary_c_root/binary_c/src -I/home/david/projects/binary_c_root/robs_extra_packages/include -I/home/david/projects/binary_c_root/binary_c/src


Building shared library for custom logging with (binary_c.h) on administrator-XPS-15-7590

Executing following command to compile the shared library:
cc -fPIC -std=gnu18 -DOPERATING_SYSTEM=linux -DLINUX -DPOSIX -DLARGEFILE_SOURCE -DMEMORY_ALLOCATION_MODEL=MEMORY_ALLOCATION_NATIVE -UALIGNSIZE -fstrict-aliasing -Wstrict-aliasing -g -Wno-sizeof-pointer-div -Wpedantic -Wshadow -Wno-variadic-macros -fstack-protector-all -rdynamic -fsignaling-nans -march=native -mtune=native -frounding-math -fno-stack-protector -ffloat-store -D__ACCURATE_BINARY_C__ -fno-finite-math-only -fasynchronous-unwind-tables -export-dynamic -O0 -DCPUFREQ=4500 -DBINARY_C_SRC=/home/david/projects/binary_c_root/binary_c/src -DPGO_off -DBINUTILS_VERSION=2.34 -DBFD_2_33 -D_FILE_OFFSET_BITS=64 -D__HAVE_LINK_H -D__HAVE__VA_OPT__ -D__HAVE_GNU_QSORT_R -D__HAVE_NATIVE_EXP10 -D__HAVE_POSIX_FADVISE -DFPU_CONTROL -D__HAVE_ATTRIBUTE___RESTRICT____ -D__HAVE_ATTRIBUTE_ALLOC_SIZE__ -D__HAVE_ATTRIBUTE_AUTO_TYPE__ -D__HAVE_ATTRIBUTE_BUILTIN_EXPECT__ -D__HAVE_ATTRIBUTE_BUILTIN_ISNAN__ -D__HAVE_ATTRIBUTE_CONST__ -D__HAVE_ATTRIBUTE_DEPRECATED__ -D__HAVE_ATTRIBUTE_GNU_PRINTF__ -D__HAVE_ATTRIBUTE_HOT__ -D__HAVE_ATTRIBUTE_MALLOC__ -D__HAVE_ATTRIBUTE_NONNULL__ -D__HAVE_ATTRIBUTE_NORETURN__ -D__HAVE_ATTRIBUTE_PACKED__ -D__HAVE_ATTRIBUTE_PURE__ -D__HAVE_ATTRIBUTE_RETURNS_NONNULL__ -D__HAVE_ATTRIBUTE_UNUSED__ -DGIT_REVISION=6895:20230518:1fb08c5af -DGIT_URL=git@gitlab-dhendriks:binary_c/binary_c.git -DGIT_BRANCH=releases/2.2.4 -D__HAVE_LIBC__ -D__HAVE_LIBCFITSIO__ -D__HAVE_LIBGSL__ -I/home/david/gsl/include -D__HAVE_LIBGSLCBLAS__ -D__HAVE_LIBDL__ -D__HAVE_LIBPTHREAD__ -D__HAVE_LIBUUID__ -D__HAVE_LIBZ__ -D__HAVE_LIBBACKTRACE__ -D__HAVE_LIBBZ2__ -D__HAVE_LIBCDICT__ -D__HAVE_LIBM__ -D__HAVE_LIBRINTERPOLATE__ -D__HAVE_BACKTRACE_H__ -D__HAVE_IEEE754_H__ -D__HAVE_DRAND48__ -D__HAVE_HSEARCH_DATA__ -D__HAVE_MALLOC_H__ -D__HAVE_SETITIMER__ -D__HAVE_HAS_INCLUDE -D__HAVE_PKG_CONFIG__ -D__HAVE_VALGRIND__ -D__SHOW_STARDATA__ -D__DIFF_STARDATA__ -D__HAVE_7Z__ -D__HAVE_BZCAT__ -D__HAVE_ZCAT__ -O0 -shared -D_SEARCH_H -L/home/david/projects/binary_c_root/binary_c/src -L/usr/lib/x86_64-linux-gnu -L/home/david/gsl/lib -L/home/david/projects/binary_c_root/robs_extra_packages/lib -L/home/david/projects/binary_c_root/binary_c/src -lc -lcfitsio -lpthread -lgsl -lgslcblas -lm -ldl -luuid -lz -lbacktrace -lbz2 -lcdict -lrinterpolate -lbinary_c -o /tmp/binary_c_python-david/notebooks/notebook_comenv/custom_logging/libcustom_logging_c2c13e8068bc4385b960647d0028913e.so /tmp/binary_c_python-david/notebooks/notebook_comenv/custom_logging/custom_logging.c -I/home/david/projects/binary_c_root/binary_c -I/home/david/projects/binary_c_root/binary_c/src -I/home/david/projects/binary_c_root/robs_extra_packages/include -I/home/david/projects/binary_c_root/binary_c/src
loading shared library for custom logging
loaded shared library for custom logging.         custom_output_function is loaded in memory at 140439734821610
Write grid code to /tmp/binary_c_python-david/notebooks/notebook_comenv/binary_c_grid_9826c6516bf44709b625ec9ff93c43ed.py [dry_run = True]
Doing a dry run of the grid.
Grid has handled 729 stars with a total probability of 0.0645564

**********************************
*             Dry run            *
*     Total starcount is 729     *
* Total probability is 0.0645564 *
**********************************

setting up the system_queue_filler now
Write grid code to /tmp/binary_c_python-david/notebooks/notebook_comenv/binary_c_grid_9826c6516bf44709b625ec9ff93c43ed.py [dry_run = False]
Signalling processes to stop

95/729  13.0% complete 18:25:37 ETA=    1.1m tpr=4.22e-01 ETF=18:26:44 mem:754.9MB M1=1.3 M2=1.1 P=3.2e+2

162/729  22.2% complete 18:25:48 ETA=    1.4m tpr=6.09e-01 ETF=18:27:14 mem:516.0MB M1=1.6 M2=1 P=2.3

287/729  39.4% complete 18:25:58 ETA=    1.1m tpr=6.09e-01 ETF=18:27:05 mem:516.6MB M1=2 M2=1.5 P=6.1e+3

426/729  58.4% complete 18:26:08 ETA=   31.9s tpr=4.21e-01 ETF=18:26:39 mem:518.6MB M1=3 M2=1.6 P=44

554/729  76.0% complete 18:26:18 ETA=   16.0s tpr=3.65e-01 ETF=18:26:34 mem:520.0MB M1=3.6 M2=3.2 P=3.2e+2

677/729  92.9% complete 18:26:28 ETA=    4.6s tpr=3.51e-01 ETF=18:26:32 mem:520.5MB M1=5.4 M2=2.7 P=17

****************************************************
*                Process 2 finished:               *
*  generator started at 2023-05-18T18:25:29.908382 *
* generator finished at 2023-05-18T18:26:34.001016 *
*             total: 1 minute and 4.09s            *
*     of which 1 minute and 3.85s with binary_c    *
*                  Ran 188 systems                 *
*       with a total probability of 0.0166562      *
*         This thread had 0 failing systems        *
*       with a total failed probability of 0       *
*   Skipped a total of 0 zero-probability systems  *
*                                                  *
****************************************************


****************************************************
*                Process 3 finished:               *
*  generator started at 2023-05-18T18:25:29.919454 *
* generator finished at 2023-05-18T18:26:34.127644 *
*             total: 1 minute and 4.21s            *
*     of which 1 minute and 3.94s with binary_c    *
*                  Ran 185 systems                 *
*       with a total probability of 0.0160873      *
*         This thread had 0 failing systems        *
*       with a total failed probability of 0       *
*   Skipped a total of 0 zero-probability systems  *
*                                                  *
****************************************************


****************************************************
*                Process 0 finished:               *
*  generator started at 2023-05-18T18:25:29.891799 *
* generator finished at 2023-05-18T18:26:35.859784 *
*             total: 1 minute and 5.97s            *
*     of which 1 minute and 5.74s with binary_c    *
*                  Ran 185 systems                 *
*       with a total probability of 0.0168338      *
*         This thread had 0 failing systems        *
*       with a total failed probability of 0       *
*   Skipped a total of 0 zero-probability systems  *
*                                                  *
****************************************************


****************************************************
*                Process 1 finished:               *
*  generator started at 2023-05-18T18:25:29.898577 *
* generator finished at 2023-05-18T18:26:36.222228 *
*             total: 1 minute and 6.32s            *
*     of which 1 minute and 6.11s with binary_c    *
*                  Ran 171 systems                 *
*       with a total probability of 0.014979       *
*         This thread had 0 failing systems        *
*       with a total failed probability of 0       *
*   Skipped a total of 0 zero-probability systems  *
*                                                  *
****************************************************


************************************************************************
*         Population-9826c6516bf44709b625ec9ff93c43ed finished!        *
*                  The total probability is 0.0645564.                 *
*  It took a total of 1 minute and 6.85s to run 729 systems on 4 cores *
*                        = 4m 27.41s of CPU time.                      *
*                     Maximum memory use 754.891 MB                    *
************************************************************************

No failed systems were found in this run.

After the run is complete, some technical report on the run is returned. I store that in analytics. As we can see below, this dictionary is like a status report of the evolution. Useful for e.g. debugging. We check this, and then set about making the plot of the orbital period distributions using Seaborn.

[7]:
print(analytics)
{'population_id': '9826c6516bf44709b625ec9ff93c43ed', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.0645563923306419, 'total_count': 729, 'start_timestamp': 1684430729.8425448, 'end_timestamp': 1684430796.6938975, 'time_elapsed': 66.85135269165039, 'total_mass_run': 3410.936346584559, 'total_probability_weighted_mass_run': 0.22609060418511753, 'zero_prob_stars_skipped': 0}
[8]:
# make a plot of the distributions
import seaborn as sns
import pandas as pd
import copy

from binarycpython.utils.functions import pad_output_distribution

pd.set_option("display.max_rows", None, "display.max_columns", None)


# set up seaborn for use in the notebook
sns.set(rc={'figure.figsize':(20,10)})
sns.set_context("notebook",
                font_scale=1.5,
                rc={"lines.linewidth":2.5})

# remove the merged objects
probability = { "merged" : 0.0, "unmerged" : 0.0}

# copy the results so we can change the copy
results = copy.deepcopy(population.population_results)

for distribution in ['post']:
    for logper in population.population_results[distribution]:
        dprob = results[distribution][logper]
        if logper < -90:
            # merged system
            probability["merged"] += dprob
            del results[distribution][logper]
        else:
            # unmerged system
            probability["unmerged"] += dprob

# pad the final distribution with zero
for distribution in population.population_results:
    pad_output_distribution(results[distribution],
                            binwidth)

# make pandas dataframe
plot_data = pd.DataFrame.from_dict(results, orient='columns')

# make the plot
p = sns.lineplot(data=plot_data)
p.set_xlabel("$\log_{10} (P_\mathrm{orb} / \mathrm{day})$")
p.set_ylabel("Number of stars")
[8]:
Text(0, 0.5, 'Number of stars')
../_images/examples_notebook_common_envelope_evolution_14_1.png

You can see that common-envelope evolution shrinks stellar orbits, just as we expect. Pre-CEE, most orbits are in the range \(10\) to \(1000\text{ }\mathrm{d}\), while after CEE the distribution peaks at about \(1\text{ }\mathrm{d}\). Some of these orbits are very short: \(\log_{10}(-2) = 0.01\text{ }\mathrm{d}\sim10\text{ }\mathrm{minutes}\). Such systems are prime candidates for exciting astrophysics: novae, type Ia supernovae and gravitational wave sources.

Things to try:

  • Extend the logging to output more data than just the orbital period.

  • What are the stellar types of the post-common envelope systems? Are they likely to undergo novae or a type-Ia supernova?

  • What are the lifetimes of the systems in close (\(<1\text{ }\mathrm{d}\)) binaries? Are they likely to merge in the life of the Universe?

  • How much mass is lost in common-envelope interactions?

  • Extend the grid to massive stars. Do you see many NS and BH compact binaries?

  • Try different \(\alpha_\mathrm{CE}\) and \(\lambda_\mathrm{CE}\) options…

  • … and perhaps increased resolution to obtain smoother curves.

  • Why do long-period systems not reach common envelope evolution?